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We introduce and study a dynamic transport model 

exhibiting Self-Organized Criticality. The novel concepts of our model are the probabilistic propa- 
gation of activity and unbiased random repartition of energy among the active site and its nearest 
neighbors. For space dimensionality d > 2 we argue that the model is related to d + 1 dimensional 
directed percolation, with time interpreted as the preferred direction. 
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Directed Percolation (DP) is one of the simplest and 
most recurrent models in statistical mechanics. Under 
very general guidelines (locality, scalar variable, etc.) 
Grassberger and Janssen [|| have proposed that a wide 
range of models would fall into the DP universality class. 
This conjecture has stood the test of time. An impressive 
parade of models under various disguises turned out to 
belong to the same DP class (see || for a review). How- 
ever, they all share one common feature: the activation 
probability p has to be defined in advance, and when it is 
properly fine-tuned, a phase transition takes place. Our 
aim here is to design a dynamical model in which critical 
directed percolation would occur via self-organization. 

Consider the following activation-transport problem. 
A conserved quantity called energy E can be stored at 
each site of a (/-dimensional square lattice with open 
boundaries. An input energy 8E < 1 is added to the 
fixed input site in the center of the system and this 
site is declared active. An active site can propagate 
activity under the following rules: 1) Random reparti- 
tion of energy among the active site and its 2d neigh- 
bors i.e. Ei — > Xi Y^jt^ 1 Ej> (i — 2d + 1) , where 

Xi = Ti/ X}?— i* r i' an d r i are uncorrelated random num- 
bers between and 1; 2) After repartition each of the 
above 2d + 1 sites becomes active with the probability 
given by its energy content. If Ei > 1 activation happens 
with certainty. There is no spontaneous activation of 
sites not connected to the current active site. The above 
process is repeated until no active sites are left in the 
system and then a new input is introduced. The above 
process strictly conserves energy except for the activation 
at the boundaries. Open boundaries allow the energy to 
escape. 

It is easy to convince oneself that the system eventually 
organizes itself into a critical stationary state where the 
average stored energy E < 1. In this state the energy 
input and output through the boundaries balance each 
other, on average. The avalanches in this critical state 



have no characteristic size except the size of the system. 
The above model is inspired by the celebrated Bak-Tang- 
Wiesenfeld (BTW) model of Self-Organized Criticality 
(SOC) {§. In fact it is a modified version of the SOC 
energy model introduced earlier ||. The new feature 
here is that the activation probability is given by the 
transported energy. 

In our model there exists an effective probability of 
activation of neighbors related to the average energy E, 
stored in the system. Each of the 2d + 1 sites participat- 
ing in a single energy redistribution has equal chance to 
become active at the next time step. The propagation of 
activity in time is reminiscent of the model of epidemic 
without immunization in d-dimensions Q , which is a par- 
ticular realization of d+l-dimensional bond directed per- 
colation with the preferred direction interpreted as time. 
The stationary state of our model corresponds to the crit- 
ical disease propagation. In Fig. |l| we show a snapshot 
of the typical fractal pattern formed by active sites at 
a given time step in the two dimensional version of our 
model. Below we shall demonstrate that despite all the 
above suggestive considerations, the aim of mapping our 
model to Directed Percolation is only partially fulfilled: 
in d > 2 indeed the activation clusters seem to fall into 
the standard DP class; for the special case d = 1 how- 
ever, though still critical, our model belongs to a new 
class, probably corresponding to DP with long range cor- 
related disorder. 

There are three quantities characterizing an avalanche 
of activity: its spatial size R, its temporal duration T, 
and its volume, equal to the total number of activations 
S. They are connected via two scaling relations: 



S = R D 
T = R z 



(1) 

(2) 



In the DP mapping these quantities are interpreted as the 
perpendicular and parallel size of cluster and its volume 
correspondingly. 
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We have measured S, T and R for the avalanches in 
our model. The results of the numerical simulation in 
1+1 dimension are shown in Fig. ||. The exponent D 
crosses over from D = 2.35±0.1, which is close to the 
DP value (2.32 @), to 1.3 ± 0.1. At the same time the 
exponent z crosses over from 1.7 ± 0.1 to a much lower 
value of 1.0 ± 0.1. It is interesting to note that S vs 
T plot has a fairly good scaling of the form S ~ T a 
throughout the whole region covered by our simulations. 
The exponent a = 1.35 ± 0.03, is different from its DP 
value a = D/z = 1.47 0. This suggests that in one 
dimension we are really dealing with a universality class 
different from 1+1-dimensional DP. 

In order to check further the relation between our 
model and DP we have performed simulations in 2+1 
dimensions. The results are shown in Fig. ||. Here we 
have measured D = 2.92 ± 0.05 and z = 1.73 ± 0.05 
in fair agreement with the DP values of 2.95 and 1.76 
f|, respectively. No crossover was observed for sys- 
tem sizes up to 120 x 120. The scaling between S and 
T is the cleanest and the measured value of exponent 
a = 1.675 ± 0.01 is in perfect agreement with its DP 
value of a = D/z — 1.675 ± 0.01 ||. We measured 
also the average stored energy E in critical state to be 
0.393 ± 0.005, for d=l, and 0.210 ± 0.005, for d = 2. 

By definition, when the avalanche is over the energy at 
every site in the system is smaller than unity. But while 
it is running the energy is not bounded and in principle 
can diverge. We expect that the larger is the avalanche 
size S the larger is the maximal energy E max (S) reached 
during this avalanche. Our numerics indicate a very slow 
increase of E max with S, for instance in 2+1 dimensions 
E m ax{S) increases from 1.18 to 1.75 while S spans 6 
decades. The observed increase in E max (s) seems to be 
even slower than logarithmic, we were not able to de- 
termine whether it would saturate or diverge slowly. In 
any case we can say that empirically events Ei > 1 are 
rare and the values are reasonably bounded. The over- 
whelming majority of the sites have simple probabilistic 
interpretation since Ei < 1. 

Another difference between our model and DP is that 
while p is a uniform constant for DP models, in our 
model Ei is space dependent. If a particular site or 
a group of sites contain more energy than the average, 
this excess amount cannot immediately diffuse to infin- 
ity. The neighborhood in question will be more prone 
to activation, which may lead to spatial and tempo- 
ral correlation in activity. Indeed this is the case for 
d = 1 . To detect spatial correlation in Ei , we have stud- 
ied the fluctuations of the convolution of the energy vari- 
able: cr(x) = ^2i =1 (Ei — E). This method is known 
to be able to detect the presence of very weak correla- 
tions in Ei Q. For uncorrelated energies we would have 
v/ (<j(x) 2 ) r~j x ' 5 . Any deviation of the scaling exponent 
from 0.5 will indicate the presence of spatial correlation 
in Ei. In the one dimensional version of our model we 



have found this exponent to be in the range between 0.33 
and 0.40, quite far from the uncorrelated value 0.5. This 
implies that a(x) performs so-called "persistent" random 
walk, being attracted to the regions it has visited before 

At first sight our model should be compared to the 
DP model with quenched disorder pd| , where the prob- 
ability of activation p(i) is given for each site and never 
changes. In our model it is true that Ei does not update 
when the activation is absent, thus in principle has a long 
memory. Activation reshuffles the energies in the neigh- 
borhood and Ei becomes "annealed" as soon as it visits 
this site or its nearest neighbors. The conserved nature 
however implies that these are not completely fresh ran- 
dom numbers, as a genuine DP model would require. For 
the model to belong to DP universality class the excess 
of energy should be able to diffuse away fast enough. 
For d — 1 as discussed above, this is not true and the 
spatial correlation of Ei is relevant. The excellent agree- 
ment of the exponents of our model in 2 dimensions with 
2 + 1-dimensional DP suggests that in high dimensions all 
these correlations are irrelevant. We conjecture that for 
d > 2 our model can effectively diffuse away any energy 
build-up and hence it belongs to the universality class of 
DP. 

On a coarse-grained level the energy is very uniformly 
distributed. This homogeneity is remarkable especially 
taking into account that the energy input in our simula- 
tions is highly inhomogeneous: the energy is injected only 
at the central site. We have checked that the exponents 
do not change when the energy is added to a randomly 
selected site each time, much in the spirit of BTW model. 
This is another manifestation of the model's ability to or- 
ganize itself to the unique critical state, independent of 
initial state or a particular type of driving. 

In traditional SOC models activation occurs via a 
threshold mechanism. Though it is possible that in 
Nature threshold activation mechanisms play an impor- 
tant role, the probabilistic activation rules are surely 
more realistic: A site containing less energy than the 
threshold can still be activated by its active neighbors, 
but this occurs less likely. Our model can be called the 
probabilistic SOC. Probabilistic activation alone does not 
make any model to fall into the DP universality class. It 
is important that we use the neutral redistribution rule, 
when each of 2c? + 1 sites involved in redistribution gets 
an unbiased random share of the conserved energy. In 
traditional SOC models, like for instance in BTW-likc 
models an active site always loses its energy to its neigh- 
bors. Therefore, during one avalanche its chances to be 
activated again are smaller than for its neighbors. We 
call such redistribution rules charitable. This leads to 
an effective repulsion of currently active sites from the 
previously activated region or outward bias in the spread 
of activity. Indeed, the dynamic exponent z ~ 1.33 in 
2-dimensional BTW model with the charitable rule of 
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dynamics is smaller than our z — 1.73. That means 
that the outward expansion of activity is faster in BTW 
model due to repulsion. It is interesting to mention that 
simple balance of the velocity of moving active front to 
the total repulsion force acting on it gives the expression 
zbtw = H which agrees very well with the results 
of numerical simulations of the BTW model. Finally a 
greedy redistribution rule where currently active sites on 
average gain energy can be envisaged. It is anticipated 
|L2| that models with charitable, neutral and greedy 
rules should belong to different universality classes. The 
detailed analysis of how the above three different redis- 
tribution rules as well as threshold or probabilistic acti- 
vation influence the universality class of SOC models will 
be relegated elsewhere. 

In summary, we have introduced a simple probabilistic 
SOC model which is a modified version of the energy 
model H . The energy is introduced to the system at the 
central site, transported around through the conservative 
repartition of energy between currently active site and its 
nearest neighbors and dissipated at the boundary. The 
novel concepts of our model are: 

1) A neutral rule of repartition of energy, where each 

site involved in it gets an unbiased random share 
of the strictly conserved sum. 

2) The probability of the site to be activated by its ac- 

tive neighbor is given by its energy. 

The system organizes itself to the critical state charac- 
terized by the power law distribution of avalanche sizes. 
The propagation of activity in time within the avalanche 
is conjectured to be equivalent to a model of spread of 
disease without immunization, which is a particular re- 
alization of bond directed percolation in d + I dimension 
with preferred direction interpreted as time. 
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FIG. 1. Snapshot of activity in one realization of the avalanche process in the system of size 400 x 400. Active sites follow a 
fractal statistics 




FIG. 2. The scaling of S and T with R in one dimension for systems of different sizes. Crossover to smaller values of 
exponents D and z happens around R = 50. The inset shows the log-log plot of S vs T in the system of size L = 3000. There 
is no crossover observed and the scaling exponent a is measured to be 1.35 ± 0.03. 
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FIG. 3. The scaling of S and T with R in two dimensions for systems of different sizes. Straight lines are fits, giving D — 2.92 
and z = 1.73. The inset shows the scaling of S with T. The measured exponent is a — 1.675 ± 0.010. 
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